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Abstract 

Any spanning tree in a loopy interaction graph can be used for communicating the effect of the 
loopy interactions by introducing messages that are passed along the edges in the spanning tree. 
This defines an exact mapping of the problem on the loopy interaction graph onto an extended 
problem on a tree interaction graph, where the thermodynamic quantities can be computed by 
a message-passing algorithm based on the Bethe equations. We propose an approximation loop 
correction algorithm for the Ising model relying on the above representation of the problem. The 
algorithm deals at the same time with the short and long loops, and can be used to obtain upper 
and lower bounds for the free energy. 
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I. INTRODUCTION 



Interesting problems are usually computationally hard and it is always useful to have 
efficient and accurate approximation algorithms. Computing local probability marginals for 
an arbitrary Gibbs measure is one of these fundamental problems in statistical physics and 
computer science. For example, having an accurate estimation of the probability marginals 
is enough to solve a constrained satisfaction problem by a decimation algorithm. 

There are, of course, some exactly solvable models that can be served as starting points 
for studying the nontrivial corrections in more interesting and complicated problems. Two 
well known examples are the mean- field (MF) solution of an infinite dimensional system 
and the Bethe solution of interacting systems that possess a tree structure [2J. In the former 
case, we are concerned with the finite dimensional corrections and in the latter, which is the 
subject of this study, we are concerned with the loop corrections. 

Given a tree interaction graph, the local probability marginals can be computed by the be- 
lief propagation (BP) algorithm, which minimizes the Bethe free energy by passing messages 
(or cavity marginals) along the graph edges (3-5]. The loopy belief propagation algorithm 
is an extension of the BP algorithm to loopy interaction graphs to find a local minimum of 
;he Bethe free energy, which could be larger or smaller than the exact free energy; see Refs. 
6-^| for some attempts to construct a convex free energy approximation and ensure the 
algorithm convergence. Nevertheless, the main strategy to deal with the loops is to group 
the variables in larger regions to eliminate some loops of given length scale as in the cluster 



variational method (9), [l0( or generalized BP 



111 ] . In the extreme limit we have the junction 



tree method [l2j, where the regions are chosen large enough to get a tree interaction graph. 
There are other algorithms that try to approximate a probability measure with simpler ones 



having tree structures [ 1 3M 1 5 ] 



The local marginals provided by the loopy BP algorithm are not necessarily consistent 
beyond the two-point correlations, and one way to improve the algorithm is to increase the 
range of consistent mare; , any fixed point of the loopy BP 

equations can be used to obtain a loop expansion of finite but exponential number of terms, 
starting from the loopy BP contribution [18| . In particular, this allows to show that for 
a class of attractive (or ferromagnetic) models the loopy BP algorithm provides an upper 



bound for the free energy 



191 ] . Finally, there are some efforts to construct field theories 



expanded around the Bethe solution [20l. |21|. 



Here we present an approximation loop correction algorithm that is based on the follow- 
ing observation: In loopy graphs, a global quantity can be computed locally by decomposing 
the computation into smaller ones distributed among different elements in the graph and 
collected by messages that are passed along the edges of any spanning tree 22, [23] . Consider 
the Ising model on an arbitrary graph and one of its spanning trees. Given a spin configu- 
ration, one can provide to each spin the effective field originated from the loopy interactions 
by passing some messages through the spanning tree. The messages are updated at each 
node to collect the effective fields coming from different parts of the tree; see Fig. [TJ In this 
way, we obtain an exact mapping of the spin configurations on the loopy graph onto the 
larger space of the spin and effective field configurations on the spanning tree. Of course, 
this mapping does not change the problem complexity but it offers some approximation loop 
correction algorithms relying on the above representation. 

In the following we discuss more about the details and write the loop correction equations 
for the Ising model. Then, we study some approximations to reduce the algorithm complexity 
by considering only a relevant subset of the loopy interactions and treating the other ones 
in a mean-field approximation. This provides an upper bound for the free energy. We also 
obtain a lower bound for the free energy by a convex combination of the loopy interactions. 



II. LOOP CORRECTIONS IN THE ISING MODEL 

Consider the Ising model with Hamiltonian H = — J2iL\ Bi s i ~ J2(ij)£E Jij s i s j 011 the 
interaction graph G = (V, E) with the set of nodes V and edges E. Given a spanning tree 
T = (V,Eq) and spin configuration s G {— 1, +1}^, we rewrite the Hamiltonian 

N 

H=-^2(Bi + BT) Si - J a s i a 3> W 

i=l (ij)eEo 

with Bj = \ Yl l j^di\a Q i^ i 3 s r Here di and d i denote the neighborhood set of i in G and 
T, respectively. The set E \ E defines the set of loopy interactions with respect to the 
spanning tree T. We are going to write the local fields Bj in terms of some messages that 
are propagated through the spanning tree. For each loopy interaction (Ik) G E \ Eq we 
introduce messages h^j, on edges (ij) in the unique path I -H- k = (I, . . . , . . . , k) 
connecting I and k on T; see Fig. [TJ In general, the messages on edge (ij) G E satisfy the 
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FIG. 1. Propagating the effect of loopy interactions through a spanning tree. Here, the loopy 
interaction Jki^k^i is replaced with the two effective fields h\_^ t = l/2JjyO"& and h k g ^ k = 1/2 Jkicri, 
which are locally determined by the messages that are passed along the edges in the spanning tree. 



following equation: 

h\^= hU i + l -Jus i {l-8 l , j ) = hU r (2) 

k£d i\j 

Consequently, we can write Bj = J2jea i ^j-H- Note that on each directed edge (i 4j)gT 
we have a vector of messages hi-yj = {h\^\l G Tj-^}. The cavity tree is denned 

recursively by i U {Tk^i\k e d$i \ j}. 

Finally, we write the partition function as 

z = j2 e ~ H = e / n j h n ^e^ (fli+E w^+s W)£ab ^ (3) 

s s i=l,...,N j&doi 

where the indicator function iff = rX/ 6 9 i ILeT^,- ^(^i^j ~~ H-^-j) ensures that the sum over 
the messages is one when the messages satisfy the equations h\_>j, otherwise it is zero. 

Now, the interaction graph is a tree and we can compute the free energy and the lo- 
cal marginals by the Bethe equations. These are self-consistent equations for the cavity 
marginals //^-(si; hy), i.e., the probability of having spin Si and messages = (hi->j, hj^i) 
in absence of node j. The equations governing the cavity marginals, called belief propagation 
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equations [3|], can be written as 



k£doi\j k£doi\j \ s k / 



• (4) 



There is, of course, one and only one solution to the BP equations that can be found by 
iteration starting from the leaves. The Bethe free energy is computed as F = Yli A-Fi — 
S(ij)eE A-Fy, where AFi and AFj are the local free energy shifts by adding node i and 
link (ij) to the interaction graph [5j, i.e., 



e- Ap > = e / n ^ x f e (Bi+E «« i,,u)8i n (j2 eJiksisk ^ s ^n > 

Si kGdoi k£do'i \ s k / 



(5) 
(6) 



Similarly, one can compute the local marginals Hi((Ji) and ^(a^aj). 



III. CONSIDERING A SUBSET OF THE LOOPY INTERACTIONS 

The time and memory complexity of the above algorithm (in the worst case) grow expo- 
nentially with the number of loopy interactions and we have to resort to some reasonable 
approximations. Here, we focus on a class of approximations that work with a subset of the 
loopy interactions and preserve the upper bound property of the free energy. The approxi- 
mation performance then depends on the structure of the spanning tree and the subset of 
the loopy interactions. 

In practice, one can start with a maximum weight W = Y^Uj)eT \ spanning tree and 
add the loopy interactions one by one according to some criterion. Then at each step one 
obtains an upper bound for the free energy of the original interacting system after adding the 
energy contribution of the discarded loopy interactions; see Fig. [2j The problem complexity 
would depend on the number of nonzero vector elements in the hi->j. Let us define X^j 
as the set of spins in the cavity tree Tj_»j that interact by loopy interactions with some 
spins in the cavity tree Tj_h. Then, the vector h^j takes 2' Xl ^ J ' values corresponding to 
the effective fields for different spin configurations in Tj_y that are relevant for the spins in 
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FIG. 2. The exact free energy F and magnetization difference Am = \ m \ ~ m i xact I ) /N vs the 
number of the loopy interactions added to the maximum spanning tree in a 2D square lattice of 
size iV = 4 x 5. We take random Gaussian fields Bi of mean zero and variance h 2 , and random 
Gaussian couplings Jij of mean zero and variance J 2 . The points are averaged over 500 independent 
realizations of the random fields and couplings. 

Let us start with a mean- field approximation of the loopy interactions, i.e., H = 
— ^2f =1 {Bi + (Bj))si — X}(r/)e£ Jij s i s j-> which results to the following BP equations: 

^•(sOoce^^W 7 *^ J] (Xy««'*/^( Sfc ) ) ■ (7) 

ke8 i\j \ s k / 

The magnetizations m, are determined self-consistently by the local marginals 

m t = L j2 s i e{Bi+ ^ ke9iXd0iJikmk)Si II (^ eJ " S!S ' c ^(^)) , ( 8 ) 

1 Si kedoi \ s k / 

where Zi is a normalization constant. 

Suppose we have partitioned the set of the loopy interactions into two subsets Cbp and 
Cmf that are to be treated exactly and in the MF approximation, respectively. For example, 
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the subset Cbp could contain all the loopy interactions within the rth neighborhood of the 
nodes in T which we denote by C r . For r = we recover the above mean- field approximation. 
Then, we rewrite the BP equations as 

k£doi\j 

x J] ^e J ^ s ^ k ^(s k ;h ik (r))Y (9) 

k£d i\j V Sfc / 

where hi->j(r) is the vector of messages going to the nodes k G Tj^i with distance da- < r. 
And Cmf^) is the set of the loopy interactions involving % that are treated in the MF 
approximation. 

Similarly, one can write the self-consistent equations for the magnetization and solve 
them by iteration starting from random initial cavity marginals and magnetizations. At the 
fixed point the cavity marginals can be used to find an upper bound for the free energy. 
There is, of course, no guarantee that the algorithm converges after introducing the MF 
approximation. 

Note that for an arbitrary subset Cbp, the computation time scales as N2 Xmax (K max X max ) 2 
with K max the maximum degree in the spanning tree, X max is the maximum |X^j|, and 
X,max — max, 

IE jedoiXj-+i\- We can indeed control the algorithm complexity by adding the 
loopy interactions as long as X max is less than a given value. 

In Figs. [3] and H] we compare the algorithm performance with the loopy BP algorithm in a 
two-dimensional (2D) square lattice. For the subset Cbp we have chosen the first Ni loopy 
interactions of largest magnitudes plus the short loopy interactions in the local neighborhood 
of radius r, such that < X max for all the directed links. 

IV. A CONVEX COMBINATION OF THE LOOPY INTERACTIONS 

One can use the above algorithm to obtain some lower bounds for the free energy. The 
free energy is a concave function of the fields and the couplings, therefore, _F[B,J] > 
P g F g \Q 9 , J 9 }, where P g is a probability measure and we need to have B = J2 g Pg^ g 



and 3 = E q P 9 J 9 



'he couplings J 9 are chosen such that the free energy F„ can be com- 



puted exactly [14|, |l5| . Let us assume that the interaction graphs G g defined by J 9 are loopy 



graphs of complexity less than X max with respect to a spanning tree T„. Finding the optimal 




FIG. 3. The free energy F and magnetization m vs the strength of the couplings in a 2D square 
lattice of size N = 50 x 50 at zero external fields Bi = with ferromagnetic couplings Jjj = J. The 
loopy BP algorithm and the exact thermodynamic solution are compared with the loop correction 
(LC) algorithm. The LC results have been obtained with a random spanning tree for different values 
of r, giving the length scale of the short loops, and X max , giving the computational complexity of 
the algorithm determined by the structure of the loopy interaction graph. 



lower bound is difficult but any consistent set of the couplings and the measure P g give a 
lower bound for the free energy. Here, we take the maximum weight spanning tree T, and 
partition the set of the loopy interactions C into subsets C r and {C g \g = 1, . . . ,J\f}, and 
set G g = T + C r + C g . Moreover, we set B 9 = B and Jf- = Jj i 5 fe)eT+£r + (Jij/P g )6(ij) e c g , 
where P g = W g / '(%2 g iW g >) with W g = Yli(ij)&c g Figure S] displays the lower bounds 

that we obtain in this way for the free energy in a 2D square lattice with random fields and 
couplings . 
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FIG. 4. The magnetization difference Am (relative to the Monte Carlo results) and the free 
energy upper and lower bounds vs the strength of the couplings in a 2D square lattice with random 
Gaussian fields Bi of mean zero and variance h 2 , and random Gaussian couplings Jij of mean zero 
and variance J 2 . The loopy BP algorithm is compared with the loop correction (LC) algorithm 
based on the maximum spanning tree and for different values of r and X max . The data points in the 
top panel are averaged over at least 10 independent realizations of the random fields and couplings. 
The lower and upper bounds in the bottom panel have been obtained for a single instance of the 
problem whereas the loopy BP results are averaged over at least 10 independent realizations. 



V. CONCLUSION 



In summary, we introduced an approximation message-passing algorithm to compute loop 
corrections in the Ising model. The approximation works with a relevant subset of the loopy 
interactions and uses the mean-field approximation to deal with the other loopy interactions. 
Obviously, the algorithm in this form is more suited to systems with strongly heterogeneous 
couplings. It would be very useful to have other approximation algorithms which treat all 
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the loopy interactions in the same manner. And finally, the algorithm can be used to obtain 
better lower bounds for the free energy than the naive one that we presented here. 
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